function [a, alpha_j, beta_j, mu_j, N_cel_j, G_array, TCP_array,  m1_array, m2_array, m3_array] = calc_TCP_metrics(dir_curves, n_ref, mu, alpha, beta, N_cel, sigma, N_stochastic, a_ref)
    %% Calculate TCPs and metrics for a set of TACs
    % dir_curves = directory to the .mat file that contains the TACs
    % n_ref = index of the TAC to take as the reference curve (default = 1)
    % mu = sub-lethal damage repair rate (1/h)
    % alpha, beta = mean values for the radiobiological parameters (1/Gy and 1/Gy^2, respectively)
    % N_cel = mean number of cells inside the tumor
    % sigma = standard deviation for parameter distribution generation (heterogeneity)
    % N_stochastic = number of random samples to generate (heterogeneity)
    % a_ref = initial choice for a (arbitrary)

    % load the fitted TACs from a .mat file
    curves_data = load(dir_curves);  
    curves = fieldnames(curves_data);
    
    % calculate the activity-dose rate scaling factor via bisection algorithm (for TCP_ref = 50%)
    TAC_ref = curves_data.(curves{n_ref});
    [a, alpha_j, beta_j, mu_j, N_cel_j] = bisection_TCP50(alpha, beta, mu, N_cel, sigma, TAC_ref, N_stochastic, a_ref);

    % initialization
    TCP_array = zeros(length(curves), 1);
    m1_array = zeros(length(curves), 1);
    m2_array = zeros(length(curves), 1);
    m3_array = zeros(length(curves), 1);
    G_array  = zeros(length(curves), N_stochastic);

    % Compute TCP and metrics for each curve
    for i = 1:length(curves)
        i
        TAC_i = curves_data.(curves{i});
        t_i = TAC_i(:,1);
        A_i = TAC_i(:,2);
        b = sum(TAC_ref(:,2)) / sum(A_i); 
        D_dot = a * A_i * b;
        G = zeros(N_stochastic, 1);
        SF_array = zeros(1, N_stochastic);
        parfor j = 1:N_stochastic
            [D , G(j)] = calc_G(mu_j(j), t_i, D_dot);
            SF_array(j) = calc_SF(D, G(j), alpha_j(j), beta_j(j));
        end
        G_array(i, :) = G;  
        TCP_tumor = exp(-N_cel_j .* SF_array);
        TCP = mean(TCP_tumor);
        TCP_array(i) = TCP;
        m1_array(i) =  metric1(TAC_ref(:,2)*a, A_i*a);
        m2_array(i) =  metric2(TAC_ref(:,2)*a, A_i*a);
        m3_array(i) =  metric3(TAC_ref(:,2)*a, A_i*a);
    end
end

function [a,alpha_j,beta_j,mu_j,N_cel_j] = bisection_TCP50(alpha, beta, mu, N_cel, sigma, TAC, N_estocastico, a_ref)
    % Calculate the activity-dose rate scaling factor (a) via a bisection algorithm
    % that aims for a TCP around 0.5 for the reference TAC
    t = TAC(:,1);
    dt = t(2)-t(1);
    R = TAC(:,2);
    
    % generate an heterogeneous distribution of radiobiological parameters
    % for the numerical (Monte Carlo) integration
    alpha_j = normrnd(alpha, sigma * alpha, 1, N_estocastico);
    beta_j = normrnd(beta, sigma * beta, 1, N_estocastico);
    mu_j = normrnd(mu, sigma * mu, 1, N_estocastico);
    N_cel_j = round(normrnd(N_cel, sigma * N_cel, 1, N_estocastico));
    
    % initialization
    a_min = 0; a_max = 100; TCP=0;
    G = zeros(N_estocastico, 1);
    parfor j = 1:N_estocastico
        [~, G(j)] = calc_G(mu_j(j), t, R * a_ref);
    end

    disp('Bisection start')
    while TCP < 0.4995 || TCP > 0.5005
    
        a = (a_min + a_max) / 2;
        R = TAC(:,2) * a;
        D = sum(R) * dt;
	    SF_array = zeros(1, N_estocastico);
    
	    for j = 1:N_estocastico
            SF_array(j) = calc_SF(D, G(j), alpha_j(j), beta_j(j));
        end
    
        TCP_tumor = exp(-N_cel_j.* SF_array);
        TCP = mean(TCP_tumor);
    
        if TCP > 0.5
            a_max = a;
        elseif TCP < 0.5
            a_min = a;
        end

    end
end

function SF = calc_SF(D, G, alpha, beta)
    % Calculate the Surviving Fraction of a tumor with radiobiological
    % parameters alpha, beta and protraction factor G after a dose D (LQ)
    log_SF = -alpha * D - beta * G * (D^2);
    SF = exp(log_SF);
end

function [D, G] = calc_G(mu, t, R)
    % calculate the protraction factor G for a given time (t), dose rate (R) curve and sublethal damage repair rate (mu)
    suma2 = 0;
    n = length(t);
    dt = t(2)-t(1);
    D = sum(R(2:n-1) * dt); 
    for i = 2:n-1
        x = 8 * log(2) / mu / dt;
        j = max(2, floor(i - x));
        suma1 = sum(R(j:i-1) .* exp(mu * (t(j:i-1) - t(i))) ) * dt;
        suma2 = suma2 + (suma1 * R(i)) * dt;
    end
    G = (2 * suma2) / (D^2);            
end

function [d1] = metric1(D_0, D_i)
    [maximo, indice_maximo] = max(D_0);
    ind = find(D_0(indice_maximo+1 : length(D_0)) >= ( maximo * 0.05 ));
    cut= indice_maximo + ind(end);
    d1 = mean( abs(D_0(1:cut) - D_i(1:cut)))/maximo;
end

function [d2] = metric2(D_0, D_i)
[maximo, indice_maximo] = max(D_0);
ind = find(D_0(indice_maximo+1 : length(D_0))< ( maximo * 0.05 ));
cut= indice_maximo + ind(1);
d2 = mean( (D_0(1:cut) - D_i(1:cut)).^2 )/maximo^2;
end

function [d3] = metric3(D_0, D_i)
[maximo, indice_maximo] = max(D_0);
ind = find(D_0(indice_maximo+1 : length(D_0))< ( maximo * 0.05 ));
cut= indice_maximo + ind(1);
d3 = max( abs(D_0(1:cut) - D_i(1:cut)))/maximo;
end